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TWO-WAY  ANALYSIS  OF  VARIANCE  FOR  WEIBULL  POPULATIONS 


1. 0 SINGLE  CLASSIFICATION  EXPERIMENTS 

During  Fiscal  Year  1978  the  work  begun  in  1977  on  the  analysis 
of  single  classification  experiments  was  brought  to  full  fruition 
with  the  completion  of  a comprehensive  (50  pages)  report  con- 
taining theory,  tables  and  worked  examples  of  both  the  analysis 
and  design  of  single  classification  experiments  with  two-parameter 
Weibull  variates.  The  report  has  now  been  accepted  for  publication  by 
the  Journal  of  Statistical  Planning  and  Inference  [1]. 

Extensive  additional  tables  of  the  numerical  values  needed 
for  analyzing  single  classification  Weibull  experiments  were 
also  computed.  Approximately  35  hours  of  CPU  time  on  SKF’s 
IBM  370/158  computer  were  used  in  generating  these  tables. 

A small  number  of  sets  have  been  reproduced  together  with  the 
report  discussed  above  as  an  SKF  publication  and  will  be 
available  to  qualified  requestors.  Copies  have  been  sent 
to  the  Rome  Air  Development  Center  for  the  use  of  the  reliability 
analysts  headquartered  at  that  site.  Copies  have  also  been  dis- 
tributed throughout  the  SKF  world-wide  organization  with  a cover 
letter  citing  the  potential  for  savings  in  bearing  testing  that 
follow  from  adoption  of  these  methods. 

The  methods  for  analyzing  single  and  crossed  classification 
experiments  developed  under  this  contract  have  been  successfully 
applied  in  other  DOD  sponsored  work  performed  at  SKF.  For  the 
Navy,  under  contract  N00019- 76-C-0168 , we  have  examined  alter- 
native plans  for  endurance  testing  three  varieties  of  silicon 
nitride  material  in  rolling  contact. 


For  the  Army,  under  Contract  DAAK70  - 77- C-00  34  we  have  compare 
the  effect  of  six  types  of  grease  on  the  endurance  performance 
of  tapered  roller  bearings. 

In  a further  effort  to  disseminate  our  results,  a presen- 
tation was  made  to  six  staff  members  of  the  Reliability  Branch 
of  Rome  Air  Development  Center  at  Rome,  N.Y.  on  March  29,  1978. 

As  a consequence  of  that  meeting  life  test  data  on  integrated 
circuits  were  commqnicated  to  us  by  RADC/RBR  for  analysis. 

The  data  were  found  to  have  been  collected  in  grouped  form  rather 
than  as  individual  failure  times  and  were  not  therefore  amenable 
to  the  analysis  developed  under  this  contract.  A more  limited 
nonparameter  analysis  was  performed  and  communicated  to  RADC  in 
the  form  of  a short  report  [2]. 

We  have  presented  a talk  based  on  the  single  classification 
analysis  entitled  "The  Comparative  Power  of  the  Likelihood  Ratio 
and  a Shape  Parameter  Ratio  Test  for  the  Equality  of  Weibull 
Scale  Parameters"  at  the  Joint  Statistical  Meetings  in  San  Diego, 
Calif.,  on  August  17,  1978.  a copy  of  the  documentation  dis- 
tributed at  the  meeting  is  given  in  Appendix  I^j 

Two  important  new  areas  of  investigation  have  been  opened 
up  during  this  contract  year,  namely,  1)  Weibull  regression  and 
2)  location  parameter  estimation  and  inference. 

2.0  WEIBULL  REGRESSION 

The  regression  problem  is  applicable  to  situations  where 
life  tests  are  conducted  at  each  of  several  levels  of  a variable 
arbitrarily  termed  a "stress".  In  applications  the  "stress" 
could  be  a voltage,  temperature,  load,  etc.  The  Weibull  scale 
parameter  or  "characteristic"  life  is  assumed  to  vary  inversely 
with  a power  of  the  stress. 
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A numerical  scheme  has  been  developed  for  the  joint  esti- 
mation by  the  method  of  maximum  likelihood  of  the  Weibull  shape 
parameter  and  the  stress-life  exponent. 


The  method  has  been  implemented  in  a computer  program 
"REGSIM'Vto  produce,  using  Monte  Carlo  methods,  the  distribution 
of  the  statistical  quantities  needed  for  setting  confidence  limits 
on  the  Weibull  shape  parameter,  the  load-life  exponent  and  a 
fixed  percentage  point  of  the  life  distribution  at  any  desired 
stress  level.  The  stress  levels  at  which  inferences  are  to  be 
made  may  be  outside  the  range  of  stresses  at  which  life  tests 
were  conducted.  The  results  are  thus  applicable  to  the  analysis 
of  "accelerated"  tests,  a common  practice  in  the  engineering  and 
physical  sciences. 


In  addition  to  REGSIM  a program  called  "REGEST"  has  been 
developed  to  perform  estimation  on  sets  of  data  taken  at  speci- 
fied stress  levels.  A paper  has  been  prepared  entitled 


"In  Terence 
and  a full 
techniques 


in  Weibull  Regression".  It  contains  theoretical  details 
numerical  illustration  of  the  application  of  the 
jA^  copy  is  included  as  Appendix  1 1 


3.0  LOCATION  PARAMETER  ESTIMATION  AND  INFERENCE 


The  Weibull  location  or  threshold  parameter  is,  in  the  life 
testing  context,  the  value  prior  to  which  failure  cannot  take 
place.  In  ordnance  applications  the  threshold  parameter  might 
represent  the  time  (after  arming)  before  which  a bomb  cannot 
explode . 


It  is  important  to  be  able  to  estimate  and  set  a lower 
bound  for  the  location  parameter  based  on  a set  of  observations. 

We  have  during  the  current  contract  year  explored  and  found 
feasible,  a technique  for  estimating  and  setting  lower  bounds  on 
the  Weibull  scale  parameter.  A program  called  l.OCEST  was  developed 
to  perform  this  analysis  on  a set  of  data.  The  program  LOCSIM 
develops,  by  Monte  Carlo  sampling,  the  tabular  data  needed  for 
implementing  the  estimation  procedure  and  optimally  selecting 
t Ifl  o ■:  If  r ■ ' • 1 i ' i V ■ i ” * ’ v c..  ' ' > ’ . 
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Several  production  runs  have  been  made  using  LOCSIM  to  try 
to  determine  a rule  for  choosing  the  early  order  statistic  at 
which  to  make  an  estimate  of  the  Weibull  slope  for  best  (most 
powerful)  detection  of  a non  zero  location  parameter.  The  fifth 
order  statistic  works  best  in  all  cases  examined. 

Below  r^  = 5 overflow  problems  were  encountered  in  the  iterative 
estimation  routines  caused  by  occasional  large  values  of  the 
estimated  shape  parameter.  We  have  incorporated  the  closed  form 
expression  for  the  ML  shape  parameter  estimate  applicable  when 
censoring  at  the  second  failure.  This  permits  trouble  free 
solution  at  ij=2  and  confirms  that  there  is  an  optimum  value  of 
r^  > 2.  Thus  the  optimum  r^  has  been  demonstrated  to  satisfy 
2 < r j < 5 in  all  cases  examined. 

Null  and  non-null  runs  have  now  been  made  in  which  the 
distribution  of  the  test  statistic  for  evaluating  the  location 
parameter  was  determined  with  r^  = 5 for  various  sample  sizes  (n) 
and  censoring  amounts  (r) , specifically:  n=10,  r=6,7 ,8, 9,10, 
n=20,  r=5, 10, 12, 15, 18, 20,  n=2S,  r=10, 14 , 18 , 20, 25,  n = 30,  r=10,15, 
20,25,30. 

An  additional  set  of  runs  was  made  with  r^  = 2 and  n=10, 
r-5,3,6,7,8,10,  n = 15,  r=5, 10,11, 12 ,15. 

This  is  sufficient  information  to  warrant  the  preparation  of 
of  a paper.  This  will  be  undertaken  early  in  the  next  contract 
year. 


A summary  of  the  software  developed  to  date  under  this 
contract  follows. 
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SOFTWARE  SUMMARY 


PROGRAM 

NAME  FUNCTION 

WEIBEST  Analyzes  k sets  of  Weibull  data.  Computes 

individual  shape  parameters,  3 percentiles 
and  printer  plot.  Calculates  shape  para- 
meter ratio  test  statistic,  likelihood  ratio 
test  statistic  and  ratio  of  largest  to 
smallest  individual  group  shape  parameter 
estimates . 

WEIBSIM  Calculates  critical  values  of  SPR  and  LR 

test  statistics  and  distributions  for 
interval  estimation  of  percentiles  and 
shape  parameter. 

REGEST  Calculates  exponent  of  power  function  rela- 

tion between  scale  parameter  and  stress, 

3 percentage  points  at  each  stress  level 
and  the  shape  parameter. 

REGSIM  Simulates  the  distribution  of  5 pivotal 

functions  needed  for  setting  single  and 
joint  confidence  intervals  in  Weibull 
regression. 

FACSIM  Calculates  the  distribution  of  functions 

required  for  testing  significance  of  row 
and  column  effects  in  factorial  experiment 
with  no  interaction. 

LOCEST  Calculates  the  test  statistic  for  location 

parameter  and  the  median  unbiased  and  lower 
951  limit  for  location  parameter. 

LOCSIM  Computes  null  distribution  of  shape  para- 

meter estimates  based  on  first  and  r 
order  statistics  in  samples  of  size  n. 
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THE  COMPARATIVE  POWER  OF  THE  LIKELIHOOD  RATIO 
AND  A SHAPE  PARAMETER  RATIO  TEST  FOR  THE 
EQUALITY  OF  WE  I BULL  SCALE  PARAMETERS 

John  I.  McCool , SKF  Industries,  Inc. 

King  of  Prussia,  PA  19406 


1.  INTRODUCTION 

We  consider  k two-parameter  Weibull  populations  having  a common 
shape  parameter  6 and  possibly  different  scale  parameters  The 

1 K 

CDF  of  the  i-th  population  is  written: 

F(x)  = 1 - exp  {-  [x/ni]P>  (1) 


Given  samples  of  size  n from  each  population,  each  type  II  censored 
at  the  r-th  failure,  the  maximum  likelihood  estimate  of  3 when  the  n- 

A 1 

are  presumed  to  differ,  is  the  solution  3i  of: 


1 ^ 

i 1 i 

k i-1  3=1 


n 3 1 n 3 1 

E xi(j)logxi(j)/_EiXi(j)] 


= 0 


(2) 


where  q ) is  the  j - th  order  statistic  in  the  sample  from  the  i-th 
population . 


The  ML  estimate  of  n-  is 

1 /V 


ni  = {1r  j^iU) 


n 3 1 
E x- 


l1/B 


(3) 


When  the  populations  have  a common  scale  parameter  n-=n  , 

A 1 O 

parameter  estimate  is  the  solution  3Q  of: 


76 


k n 3o 

J *iCj)l08*i(j) 

o * k £ E l0*‘i[i]  - — J~ 

° i=l  j = l 1UJ  k n 


k r 


= 0 


1 E xih) 


The  ML  estimate  of  nQ  is 


k n 3, 


Mi  A x‘«>/rk) 


~)1/eo 


the  ML  shape 


(4) 


(5) 
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TESTING  EQUALITY  OF  SCALE  PARAMETERS 

We  consider  two  tests  of  the  hypothesis  H :n  =n  a shape  para- 
meter  ratio  (SPR)  test  based  on  the  ratio  81/80  of  the  two  shape  parameter 
estimates  of  Eqs . (2)  and  (4)  and  a test  ba^ed  on  the  likelihood  ratio 
X using  as  test  statistic  the  monotonically  transformed  value  -21ogX 
expressible  as, 

- 21ogX  = 2 rklog(B  1/80)  + 2 E E [log(-i^  ) 8 1 - logC-i^  ) 00  ] (6) 

i=l  j=l  hi  n0 

The  distributions  of  both  test  statistics  depend  only  upon  n,  r and  k 
under  the  null  hypothesis. 

A A 

Table  1 gives  90  and  95%  points  of  81/80  and  -21ogX  determined  by 
Monte  Carlo  sampling  for  various  n,  r and  k with  n ranging  from  5 to  30 
and  k,  from  2 to  10. 

The  degree  to  which  -21ogX  approaches  the  asymptotic  x2  distribution 
based  on  k-1  degrees  of  freedom  may  be  judged  by  their  respective  95-th 
percentiles,  listed  in  columns  7 and  8 of  Table  1. 

POWER  OF  HYPOTHESIS  TESTS 

The  power  of  the  SPR  test  against  various  alternatives  has  been 
found  to  depend  only  upon  t.he  value  of  the  symmetric  function  of  the  q.  : 

<h  = £ lognf  [n?/  En?  - X/k]  (7) 

i=l  1 1 i=l1 

<+>  1 is  scale  invariant.  When  the  ' s are  equal,  <J>i  = 0,  and  when  the 
are  not  all  equal,  4>i>0. 

The  power  of  the  LR  test  depends  on  <J>  1 anc  another  symmetric, 
scale  invariant,  postive  function  of  the  n^'s  given  by 

k k 

<J>2  = £ log  [Enf/knf ] (8) 

i=l  i-11  1 

For  equal  $ 1 the  power  of  the  LR  test  may  be  slightly  inferior  or 
slightly  superior  to  the  SPR  test,  depending  upon  the  value  of  <{>2.  In 
general,  SPR  is  best  against  a single  aberrant  i.e.,  when  all  n^'s  are 
equal  but  one,  and  k>2;  the  LR  test  is  superior  against  more  "diffuse" 
alternatives  in  which  several  ' s differ. 
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For  fixed  r the  power  of  both  tests,  against  the  single  aberrant 
alternative,  decreases  with  n and  increases  with  k. 

Figure  1 shows,  for  k=2  and  n=30  the  value  of  the  acceptance  pro- 
bability plotted  for  various  r against  n?  (n2=1.0).  The  SPR  and  LR 
tests  are  indistinguishable  to  graphical  accuracy.  These  values  may 
be  used,  for  each  r as  a conservative  approximation  for  the  acceptance 
probabilities  against  the  single  aberrant  alternative,  when  k>2  and 
n<30. 

4 . INTERVAL  ESTIMATION  OF  THE  SHAPE  PARAMETER 

The  simulations  that  were  used  to  find  the  critical  values  of  the 
SPR  and  LR  tests  also  yielded  the  distributions  of  the  function 

A A 

v = gi/g  (9) 

This  distribution  depends  only  on  n,  r and  k.  A 100(l-a)%  interval 
estimate  of  g is  obtained  by  inverting  the  inequalities  in  the  proba- 
bility statement 

Pr  [vo/2  < 6,/S  < Vo/2]  - 1 - „ ,10) 

For  90%  confidence  the  interval  is 

®'/v0.95  < 8 < 8 1 ^ v0 , 0 5 <“> 

The  ratio  R of  the  upper  to  lower  end  of  this  interval  is  given  by 

R = V0 . 9 5 ^ V0 . 0 5 *12) 

The  quantity  R may  be  used  to  characterize  the  precision  within 
which  g has  been  determined  by  the  k censored  samples  of  size  n. 

Figure  1 shows  some  plots  of  R against  censoring  number  for  n=5,  15 
and  30  and  with  k=l  and  k=5.  For  a given  value  of  r the  precision 
worsens  with  n.  For  fixed  n precision  improves  with  r.  Considerable 
precision  is  gained,  particularly  at  low  censoring  numbers,  if  samples 
can  be  pooled  for  estimating  g,  i.e.,  if  k>  1 . 
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5 . INTERVAL  ESTIMATION  OF  A PERCENTILE  OF  THE  i-th  POPULATION 

The  100  p-th  percentile  of  the  i-th  population  is  expressible  in 
terms  of  and  8 as 


x = [-log(l-p) ]1/3*  n- 
pi  1 

Its  ML  estimate  is  correspondingly 

xp  = [-log(l-p) ]1/01  • ni 


(13) 


(14) 


The  function 

a.  x 

u s 3 1 log(  pi/x  ) (15) 

Fi 

follows  a distribution  that  depends  on  n,  r and  k only.  In  terms  of 
the  5-th  and  95-th  percentiles  of  u,  a 90%  confidence  interval  for  xp 
is  calculable  as 


exp  [■ 


U0.95/B] 


< x 


< x 


Pi 


CXP  [~uo.  05/g1 


(16) 


The  ratio  of  the  upper  to  lower  ends  of  this  confidence  interval 
is  a random  variable  through  its  dependence  on  8.  The  8-th  power  of 
the  median  value  of  this  random  variable  denoted  may  be  expressed  as 


R 


8 

0.50 


exp[(-u005  + u0>g5)/v0S0] 


(17) 


and  may  be  used  as  a measure  of  the  precision  with  which 
determined. 


x has  been 
pi 


O 

Figure  2 shows  R^  j-q  for  estimating  the  tenth  percentile  x^ 
of  any  population,  as  a function  of  censoring  number  for  some  n and  k 
values.  Again,  pooling  of  samples  believed  to  have  a common  shape 
parameter  greatly  increases  the  precision  in  the  determination  of  Xq 
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A A 

ll/ft)  -21°g>  >2(lc-l) 


n 

r 

k 

0.90 

0.95 

0.90 

0.95 

0.95 

5 

3 

2 

1.  799 

2.144  ‘ 

5. 357 

7.  235 

3.  84 

5 

3 

3 

1.788 

2.078 

7.913 

10.24 

5.99 

5 

3 

4 

1.  769 

2.009 

10.29 

12.80 

7.  81 

5 

3 

5 

1.744 

1.925 

12. 34 

15.07 

9.49 

5 

. 3 

10 

1.644 

1.747 

21.  78 

24.88 

16.9 

5 

5 

2 

1.259 

1. 366 

3.  794 

5.290 

3.  84 

5 

5 

3 

1.273 

1. 361 

5.984 

7.  804 

5.99 

5 

5 

4 

1.2  75 

1. 35  3 

8.0  76 

10.08 

7.  81 

5 

5 

5 

1.267 

1. 339 

9.843 

12.10 

9.49 

5 

5 

10 

1.246 

1.284 

17.94 

20.72 

16.9 

10 

5 

2 

1. 36  3 

1.508 

4.142 

5.782 

3.  84 

10 

5 

3 

1. 372 

1.491 

6.443 

8.  372 

5.99 

10 

5 

4 

1.36  7 

1.463 

8.467 

10.51 

7. 81 

10 

5 

5 

1.  365 

1.438 

10.39 

12.45 

9.49 

10 

5 

10 

1.328 

1. 378 

18.84 

21.58 

16.9 

10 

10 

2 

1 .105 

1.150 

3.165 

4.472 

3.84 

10 

10 

3 

1.116 

1.149 

5. 317 

6. 834 

5.99 

10 

10 

4 

1.115 

1.142 

6.981 

8.  706 

7.81 

10 

10 

5 

1.114 

1 . 1 39 

8.755 

10.55 

9.49 

10 

10 

10 

1.105 

1.121 

16.16 

18.61 

16.9 

15 

5 

2 

1.  3 74 

1.5  35 

4.090 

5.  770 

3.  84 

15 

5 

3 

1.392 

1.522 

6.393 

8.428 

5.99 

15 

5 

4 

1.392 

1.491 

8.5  79 

10.63 

7. 81 

15 

5 

5 

1.377 

1.462 

10.23 

12.44 

9.49 

15 

10 

2 

1.131 

1.190 

3.260 

4.653 

3.  84 

15 

10 

3 

1.144 

1.185 

5.  370 

6.846 

5.99 

15 

10 

4 

1.145 

1.181 

7.224 

8.962 

7.  81 

15 

10 

5 

1.143 

1.176 

8.881 

10.87 

9.49 

15 

15 

2 

1.066 

1.096 

3.101 

4.4  36 

3.  84 

15 

15 

3 

1.072 

1.095 

5.071 

6.562 

5.99 

15 

15 

4 

1.0  74 

1.091 

6.843 

8.497 

7. 81 

15 

15 

5 

1.0  74  . 

1.090 

8.477 

10. 38 

9.49 

20 

5 

2 

1.382 

1.546 

4.110 

5.  749 

3.  84 

20 

5 

3 

1.408 

1.542 

6.455 

8.469 

5.99 

20 

5 

4 

1.39  8 

1.495 

8.487 

10.51 

7.  81 

20 

5 

5 

1.389 

1.478 

10. 34 

12.41 

9.49 

20 

10 

2 

1.139 

1.197 

3.220 

4.587 

3.84 

20 

10 

3 

1.155 

1.202 

5.45  3 

7.004 

5.99 

20 

10 

4 

1.155 

1.195 

7.233 

9.120 

7.  81 

20 

10 

5 

1.155 

1.190 

9.068 

11.01 

9.49 

20 

15 

2 

1.079 

1.112 

3.088 

4.297 

3.  84 

20 

15 

3 

1.087 

1.115 

5.102 

6.583 

5.99 

20 

15 

4 

1.088 

1.111 

6.862 

8.591 

7.  81 

20 

15 

5 

1.089 

1.108 

8.603 

10.56 

9.49 

Table 

1.  90-th 

5 95-th 

Pe  rccntiles 

of 

1 


S I’  R an d I . U To s t Statistics 


A . 


- 2 log  A 


X-(k-l) 


n 

jr 

k 

0.90 

0.95 

0.90 

0.95 

0.95 

20 

20 

2 

1.048 

1.068 

3.012 

4.194 

3.  84 

20 

20 

3 

1.054 

1.069 

4.986 

6.456 

5.99 

20 

20 

4 

1.054 

1.068 

6.740 

8.460 

7.  81 

20 

20 

5 

1.053 

1.065 

8.428 

10.23 

9.49 

25 

5 

2 

1.383 

1.545 

4.003 

5.6  34 

3.  84 

25 

5 

3 

1.400 

1.522 

6.365 

8.142 

5.99 

25 

5 

4 

1.402 

1.498 

8.518 

10.40 

7.  81 

25 

5 

5 

1.  390 

1.481 

10.29 

12.42 

9.49 

25 

10 

2 

1.143 

1.199 

3.211 

4.485 

3.  84 

25 

10 

3 

1.160 

1.208 

5.403 

7.038 

5.99 

25 

10 

4 

1.161 

1.202 

7.274 

9.074 

7. 81 

25 

10 

5 

1.158 

1 .191 

8.872 

10.83 

9.49 

25 

15 

2 

1.084 

1.122 

3.  079 

4.465 

3.  84 

25 

15 

3 

1.094 

1.123 

5.116 

6.633 

5.99 

25 

15 

4 

1.096 

1.123 

* 6.992 

8.  885 

7.  81 

25 

15 

■ 5 

1.095 

1.117 

8.544 

10.52 

9.49 

25 

20 

2 

1.0  56 

1.078 

2.994 

4.250 

3.  84 

25 

20 

3 

1.063 

1.081 

4.98S 

6.453 

5.99 

25 

20 

4 

1.06  3 

1 .080 

6.  783 

8.448 

7.81 

25 

20 

5 

1.062 

1.077 

8.212 

10.15 

9.49 

25 

25 

2 

1.037 

1.054 

2.945 

4.210 

3.  84 

25 

25 

3 

1.042 

1.054 

4.894 

6.450 

5.99 

25 

25 

4 

1.042 

1.053 

6.616 

8.  330 

7.  81 

25 

25 

5 

1.041 

1.051 

8.122 

9.919 

9.49 

30 

5 

2 

1.394 

1.571 

4.111 

5.852 

3.  84 

30 

5 

3 

1.414 

1.544 

6.489 

8. 396 

5.99 

30 

5 

4 

1.411 

1.529 

8.601 

10.83 

7.  81 

30 

5 

5 

1 . 39  9 

1.492 

10.3  7 

12.67 

9.49 

30 

10 

2 

1 . 1 4 S 

1 . 212 

3.  259 

4.  704 

3.  84 

30 

10 

3 

1.  165 

1 .215 

5.452 

7.200 

5.99 

30 

10 

4 

1.16  7 

1 . 206 

7.402 

9.171 

7.  81 

30 

10 

5 

1.165 

1.199 

9.071 

11.03 

9.49 

30 

15 

2 

1.087 

1.124 

3.033 

4. 317 

3.  84 

30 

15 

3 

1.096 

1.125 

5.046 

6.549 

5.99 

30 

15 

4 

1.099 

1.12  5 

6.96  7 

8.663 

7.81 

50 

15 

5 

1.098 

1.121 

8.588 

10.51 

9.49 

30 

20 

2 

1.059 

1.083 

2.921 

4.167 

3.  84 

30 

20 

3 

1.065 

1.084 

4.844 

6.256 

5.99 

30 

20 

4 

1.067 

1.084 

6.687 

8. 405 

7.81 

30 

20 

5 

1.06  7 

1.082 

8.  361 

10.17 

9.49 

30 

25 

2 

1.044 

1 . 061 

2.9  56 

4.116 

3.84 

30 

25 

3 

1.047 

1.061 

4.753 

6.168 

5.99 

30 

2 5 

4 

1.04S 

1.060 

6.562 

8.082 

7.  81 

30 

25 

5 

1.048 

1.058 

8.128 

9.886 

9.49 

Table  1 (Continued) 
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-21ogX 

y2(k-l) 

n 

r 

k 

0.90 

0.95 

0.90 

0.95 

0.95 

30 

30 

2 

1.0  31 

1.044 

2. 890 

4.025 

3.  84 

30 

30 

3 

1.0  35 

1.044 

4.765 

6.192 

5.99 

30 

30 

4 

1.035 

1.044 

6.511 

8.123 

7.  81 

30 

30 

5 

1.0  35 

1.042 

8.121 

9.938 

9.49 

Table  1 (Continued) 
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INFERENCE  IN  WEI  BULL  REGRESSION 


ABSTRACT 

We  consider  an  experimental  situation  in  which  a response 
variable  follows  a two -parameter  Weibull  distribution  having  a 
scale  parameter  that  varies  inversely  with  a power  of  a deter- 
ministic, externally  controlled,  variable  generically  termed  a 
"stress".  The  shape  parameter  is  presumed  to  be  invariant  with 
stress.  Equations  are  formed,  from  the  results  of  type  II  cen- 
sored life  tests  conducted  at  each  of  several  stresses,  whose 
solution  yield  the  maximum  likelihood  estimates  of  the  common 
shape  parameter,  the  stress-life  exponent,  and  a general  percen- 
tile of  the  life  distribution  applicable  at  an  arbitrary  stress 
level.  A numerical  scheme  for  solving  the  equations  is  given. 

Pivotal  functions  of  the  ML  estimates  are  found  whereby 
interval  and  median  unbiased  point  estimates  may  be  calculated 
once  the  distribution  of  the  pivotal  functions  is  found  by 
Monte  Carlo  sampling.  A numerical  example  of  the  calculation  of 
point  and  interval  estimates  is  given.  The  precision  with  which 
the  shape  parameter  is  estimated  by  testing  at  several  stresses 
is  comparable  to  the  precision  applicable  to  a single  test  of 
the  same  total  sample  size.  The  precision  in  estimating  percen- 
tiles is  a maximum  near  the  middle  of  the  stress  range  at  which 
testing  was  performed  and  is,  at  that  point,  comparable  to  the 
precision  obtained  in  a single  stress  test  of  the  same  total 
sample  size. 


INFERENCE  IN  WE I BULL  REGRESSION 


1 . INTRODUCTION  AND  SUMMARY 

The  conduct  of  life  tests  at  more  than  one  level  of  an  environ- 
mental factor  known  to  affect  the  parameters  of  the  life  distri- 
bution is  a common  practice.  Ordinarily,  the  environmental  fac- 
tor is  set  at  higher  levels  of  intensity  than  the  test  item  will 
meet  in  service,  in  order  to  shorten  the  expected  time  to  failure. 
The  life  test  results  obtained  at  these  levels  are  then  extra- 
polated to  more  usual  levels  of  the  environmental  factor  by 
fitting  constants  to  a theoretical  or  empirical  relation  between 
the  factor  levels  and  one  or  more  parameters  of  the  life  distri- 
bution. (It  is  usual  to  assume  the  form  of  the  life  distribu- 
tion is  not  altered  by  the  level  of  the  environmental  factor.) 

Life  testing  having  these  aims,  lias  come  to  be  called  accelerated 
life  testing  and  an  extensive  literature  on  the  subject  has 
developed,  primarily  in  connection  with  the  testing  of  electronic 
components.  The  environmental  or  accelerating  factor  is 
customarily  called  stress,  but  may  actually  be  voltage,  load, 
temperature,  etc. 

The  work  described  in  this  paper  is  applicable  to  the  analysis 
of  censored  accelerated  life  tests  in  which  i)  a two  parameter 
Weibull  distribution  governs  at  each  stress  and  ii)  the  Weibull 
scale  parameter  varies  inversely  with  a power  of  the  stress  while 
iii)  the  shape  parameter  is  invariant  with  stress. 

The  primary  aims  of  the  present  work  are  to  find  point  and 
interval  estimators  of  i)  a quantile  of  the  life  distribution 
at  a specified  stress,  ii)  the  Weibull  shape  parameter  and 
iii)  the  exponent  of  the  stress-life  model  based  upon  the 
method  of  maximum  likelihood.  These  aims  arc  analogous  to  the 
traditional  aims  of  regression  analysis. 
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Singpurwalla  [1]  has  considered  the  case  where  the  lives  follow 
the  single  parameter  exponential  distribution  with  the  scale 
parameter  varying  inversely  with  a power  of  the  stress  and 
where  the  data  at  each  stress  level  is  censored.  He  derives 
the  maximum  likelihood  (ML)  estimates  of  the  stress-life 
exponent  and  the  constant  of  proportionality  of  the  stress- 
life  law.  He  uses  asymptotic  theory  for  setting  confidence 
limits  on  these  parameters. 

Nelson  [2]  treats  the  two  parameter  Weibull  case  but  estimates 
the  parameters  separately  at  each  stress.  He  then  estimates 
the  common  shape  parameter  as  a weighted  combination  of  the 
individual  shape  parameter  estimates  and  the  stress-life 
exponent  and  proportionality  constant  by  weighted  least  squares 
using  the  logarithmically  transformed  stress  values  as  the 
independent  variable. 

In  [3]  Singpurwalla  and  Al-Khayyal  under  the  same  assumptions  as 
consider  the  direct  use  of  the  method  of  maximum  likelihood 
to  estimate  the  common  shape  parameter  and  the  stress-life 
exponent  and  constant  and  find  the  asymptotic  covariance 
matrix  of  the  estimators  for  the  uncensored  case. 


In  a numerical  example  they  calculate  these  estimates  by  direct 
maximization  of  the  likelihood  function  using  SUMT. 


In  Section  2 of  the 
are  derived  for  the 
the  stress-life  law 
bution  at  arbitrary 


present  paper,  the  ML  estimation  equations 
joint  estimation  of  the  two  parameters  of 
and  a general  quantile  of  the  life  distri- 
s tress  S for  the  case  where  n^  items  are 
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. . ....  1 ri — 


tested  at  each  of  k stress  levels  S^ . . . S^  and  testing  continues 
to  the  occurrence  of  the  r^-th  failure  at  each  stress. 

In  Section  3,  the  ML  equations  derived  in  Section  2 are  tested 
to  determine  whether  they  reduce  to  the  ordinary  equations  when 
a single  stress  level  is  employed  and  whether  the  estimates  are 
properly  invariant  when  the  units  of  the  stress  scale  are 
changed. 

Section  4 contains  the  derivation  of  pivotal  functions  which, 
given  their  distributions  for  specific  sample  sizes,  will  permit 
bias  correction  of,  and  the  setting  of  confidence  limits  on, 

1)  a general  quantile  of  the  life  distribution  at  stress  S, 

2)  the  Weibull  shape  parameter  and  3)  the  stress-life  exponent. 

Section  5 contains  the  description  of  a scheme  found  to  be 
effective  for  the  numerical  solution  of  the  likelihood  equations 
and  briefly  describes  a computer  program  for  the  Monte  Carlo 
calculation  of  the  distribution  of  the  pivotal  functions 
derived  in  Section  4. 

Section  6 is  a numerical  example  of  the  analysis  of  a set  of 
rolling  contact  fatigue  data  obtained  at  four  levels  of  the 
contact  stress. 

2.  PROBLEM  FORMULATION  AND  DERIVATION  OE  ESTIMATORS 

We  consider  a series  of  life  tests  in  which  a sample  is  tested 
at  each  of  k stress  levels.  We  suppose  that  for  i = l...k, 
n^  specimens  are  subjected  to  stress  S^  and  run  until  the 
occurrence  of  the  first  r^  failures. 


It  is  assumed  that  under  stress  the  time-to-failure  of  the 
test  item  is  a two-parameter  Weibull  random  variable  having  a 
scale  parameter  r\ ^ and  a stress -independent  shape  parameter  3. 
That  is,  the  cumulative  distribution  function  of  the  life  at 
stress  is  expressible  as 

Prob[life  <x|  S=S^]  = 1 - exp [ - (x/n^) (1) 


The  scale  parameter  is  assumed  to  vary  inversely  with  the 
Y-th  power  of  the  stress  S^,  i.e.  as: 


n . = 
1 


(2) 


where  nQ  is  a constant  representing  the  scale  parameter  at  unity 
stress . 


We  denote  by 

wUh  xi(n) 


xi(j)  the 
Xi(ri+l) 


j-th  ordered  life  achieved  at 
Xi(ni)’ 


stress 


Given  the  results  of  such  a life  test,  we  wish  to  form  the 
maximum  likelihood  (ML)  estimates  of  the  parameters  nQ,  6 and 
Y. 


The  logarithm  of  the  likelihood  function  written  for  all  k 
samples  is  found  to  have  the  form 


k k ri 

logL  = logc  + log  3 E r.  + (3-1)  E E log  x.,,,.  (3) 

i=l  1 i = l j = l lUJ 

k k ni 

-8  rjiog  nj  - 
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Equating  to  zero  the  derivative  of  logL  with  respect  to  nQ 
gives  the  following  equation  which  must  be  satisfied  by  the  ML 

A A A 

estimates  q , g and  y. 


= 'B/rio  z ri  + *n0-p-1  z z (-iai)B=  0 

11  i=l  j=l  S."Y 


-6-1 


2 2 /i(i)^  = 


The  solution  for  n in  terms  of  y and  6 is 

o 

% = { z z1  ^ ri}  1/P 

i=l  j=l  Si  Y i=l 


When  n^  = n and  r.  = r for  all  i this  equation  specializes  to 


- i k n x;m  q 1/6 

"o  ■ E 1 > 1 

° rx  i= 1 j=l  Si  Y 


Differentiating  Eq.  (3)  with  respect  to  y and  substituting 
from  Eq.  (5)  gives: 

k A ri  A v8  ioeei  v'1"'1' 

**  ■ ■ f i»- -'r,;  i 0 

if!  Si  X Xi ( j ) 

For  r^  = r,  n^  = n,  Eq . (7)  reduces  to: 


9 logL 
3 Y 


Z SyB  logS  Z [xim] 
i=  1 1 1 i = l UJ 


I loESi  - — - V S 

1 = 1 y C Y6  y x 

i=  1 1 J-*! 


■ 0 (8) 


‘ 


Finally,  differentiating  Eq . (3)  with  respect  to  3 and  using 
Eqs.  (5)  and  (7)  leads,  after  considerable  simplification,  to 
the  expression: 


r • k ^ ^ a 

1 ifl  j.El  l0SXitj)  ih  Si1,6jf1  Xi(j)  l0gXi(j) 

^ + i-r * V zr. »* 

0 E r E SYB  1 X3 

w 1 Li  1 i=i  iCj) 


= 0 


(9) 


5-  C Y3  y 3 
i-1  1 j-1  Xi«)  Sxi(5) 


For  r^  = r and  n^  = n,  Eq.  (9)  specializes  to: 

k 
r 

i i=r  1 

n 

s.YP 

t 

i=l 


k 

i i . r £ logx . > 

~ + 7F  i-l  j-1  l(j3 


z s.yB  Z x?,.. 

1 j-1 


= 0 (10) 


Following  the  simultaneous^ solution  of  Eqs.  (7)  and  (9)  for 
0 and  Y one  nay  evaluate  nQ  from  Eq.  (5).  The  ML  estimate  of 
0.  is  obtained  using  Eq.  (2)  and  the  fact  that  ML  estimates  of 
functions  of  parameters  are  just  those  functions  of  the  ML 
estimates,  i.e. 


A 


Si 


-Y 


(ID 


The  p-th  quantile  of  the  life  distribution  under  stress  may 
be  estimated  as: 

, 1/3  * 

*p.  = [i°s(]4p)]  ‘ ni  (12) 

3 . CHECKING  THE  VALIDITY  OF  THE  EQUATIONS 


As  a check  on  the  reasonableness  of  the  likelihood  equations 
given  above  we  consider  the  special  case  in  which  all  testing 
is  at  a single  stress,  i.e.  = S (all  i).  In  this  case 


the  stress-life  exponent  y should  be  unestimable  and  the  equa- 
tions for  estimating  g and  should  reduce  to  the  expression 
known  to  apply  when  a single  sample  is  considered. 


From  Eq.  (8)  with  = S one  has: 

k n R 
k SYB  logs  E Z 

i i r-  i=l  j=l  ^ J _ 

klogS ^ g - ^ J 


SYP  Z Z xB 
i=l  j = l 


(13) 


which  reduces  to: 

k log  S - k log  S = 0 

A 

i.e.,  the  equation  is  satisfied  for  any  value  of  y. 


Eq.  (10)  with  Si  = S becomes 


1 . 1 


7E  . 


k r 
l l logx . ... 


"„k  n 2 
; Y 6 ^ v 3 


8 *“  i = l i = i vi  k " b 

B 1 1 j 1 SYP  E E xr,.. 

i-1  j-1  1(J) 


E E x.,..  1 o gx  • , . 
i-i  i-i  lU)  >0)  „ „ 

— t: — — * “ u 


(14) 


The  estimate  of  nQ  becomes,  from  Eq . (6): 


* V n ~ 

- SY  f1  E E xB  }1/B 

_ S % jfi  Xi(j)} 


(15) 


The  estimate  of  , the  Weibull  scale  parameter  at  stress  S using 
Eq.  (11)  then  becomes: 


: _ {1  j ? x e 

i£-l  ill  UJ)’ 


(16) 
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Eqs.  (14)  and  (16)  arc  seen  to  be  equivalent  to  the  ML  equations 
for  estimating  n and  3 from  a single  sample  of  size  kn  having 
rk  failures  (cf.  Cohen  [4]). 


Another  check  for  reasonableness  of  the  equations  can  be  made 
by  virtue  of  the  fact  that  the  estimates  of  6,  Y and  n.  should 
be  invariant  with  respect  to  a change  in  scale  of  the  stress 
e.g.  the  estimates  should  not  be  affected  if  the  units  in 
which  stress  is  measured  are  changed  from  English  to  metric. 

To  test  whether  this  is  true,  introduce  a scale  change 


ai  = cSi  (17) 

where  c is  a constant  and  is  a multiplicatively  transformed 
value  of  S^. 


Substituting 


k 

£ logo.  - klogc 
i=l  1 


°i/c  into  Eq.  (8)  gives: 


n 


k c,  ye  .. 

k E (-i)  [logo,  - logc]  I x ni 

i=l  j = l 1 u J 

— ’ — * 


i = l 


<cT> 


i-vYe  “ .6 


o (18) 


XT 


j=l 


i(j) 


which  becomes  : 


k 

E logo.  - k logc 
i = l 1 


k E o • 
i = l 1 
k 


V3 


logo. 


n 


E o E X6 

E °i  ^ xuj) 


+ k logc  = 0 


i = l 


(19) 


Substituting  S^  = o^/ c into  Eq . (10)  gives: 


8 


k r 


Y6  " 8 


1.1  i i , tl/c  W1  4ixi(j)lo8xi(j) 

I “‘uj.-  „ '*vi";vv  ? _s '0C20) 


S **WJ-l  ‘'"-.u/c’VV*  1 4) 

i=l  1 j=l  1UJ 

Eqs.  (19)  and  (20)  are  identical  to  Eqs.  (8)  and  (10)  respectively 

except  that  replaces  S-.  Thus  the  simultaneous  solution  of 
i l ~ 

Eqs.  (8)  and  (10)  for  ? and  3 is  invariant  with  respect  to  a 
scale  change  in  the  stresses. 

A 

The  solution  for  using  = o ./c  becomes  from  Eq . (6): 


c-i,i  i 5 (WS,i/* 

rK  i=l  i=l  o i 

J i 


Eq.  (11)  then  becomes: 


. 11/  C-V(1  \ s 

1 c y rk  i=l  j=l  O.'Y 

J l 


Again,  Eq.  (22)  is  identical  to  Eq.  (11)  with  replacing  , 
so  that  the  ML  estimate  of  n-  is  unaltered  by  a scale  change 
in  stress. 


4.  PIVOTAL  FUNCTIONS 


It  is  possible  to  draw  inferences,  i.e.  set  confidence  inter- 
vals and  make  hypothesis  tests  for  a parameter  if  one  can 
determine  a function  of  the  parameter  and  its  estimate  that 
follows  a distribution  that  depends  on  sample  sice  but  not 
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upon  that  parameter  or  any  other  parameters  of  the  distribution. 
Such  functions  are  designated  pivotal  functions  and  have  been 
found  for  the  shape  parameter  and  a general  quantile  in  the 
single  sample  case  of  the  two  parameter  Weibull  distribution 
[cf.  McCool  (5)]. 

In  searching  for  pivotal  functions  for  the  present  problem  we 
use  the  strategy  employed  in  [5]  of  expressing  a Weibull  random 
variable  in  terms  of  its  population  parameters  and  a rectangular 
variable.  When  this  is  done  the  order  statistic  transforms 

to 


Xi(j)  ‘ ni  (-lo2uij  1 


where  the  variables  follow  beta  distributions  with  parameter? 
that  depend  on  sample  size  but  not  upon  or  B . 

For  simplicity  we  hereafter  restrict  attention  to  the  special 
case  of  constant  sample  size  in  which  n^  = n and  r^  = r.  The 
results  apply  to  the  unequal  sample  size  case  as  well. 


Substituting  Eq.  (23)  into  Eq . (8)  and  using  (2)  results,  after 
some  simplification,  in  the  equation: 


l logS.  - 


k ^ ^ ^ n 

k E S.e(Y"Y)logS.  £ {-logu. 
i=l  1 1 j=l  1J 

Z S,PtY'Y)  Z {-logu. . } 6/p 


r. 


Substituting  Eq.  (23)  into  Eq.  (10)  and  using  Eq.  (24)  leads  to: 

k n 


k r 


E E s?(Y'Y)<-logu ii)6/Sloglog(l/uii) 
i=l  1 = 1 1 1J  1J 


S ■ 13  J E SS(Y'Y)  [-logu..]S/B 


0 (25) 


i-lj-1 


If  one  now  writes  g(y-y)  as  (3/3)(y-y)g  it  is  seen  that  for  a 
given  set  of  u- . Eqs.  (24)  and  (25)  can  be  solved  simultaneously 

lj  A /N 

for  the  quantities  q2  g/g  and  SE(y-y)g  . In  repeated  sampling 
i.e.  different  sets  of  u^j , these  quantities  will  vary  in  a 
manner  that  depends  only  on  k,  r and  n. 


If  for  specified  values  of  k,  r and  n repeated  Monte  Carlo 
samples  were  drawn  from  a two  parameter  Keibull  population 
having  say  g = 1.0  and  y = 0,  one  could  empirically  determine 

A 

as  closely  as  desired  the  distribution  of  3 and  hence  of  q and 
the  distribution  of  y and  hence  of  s. 

Denoting  the  100  a-th  percentage  point  of  the  distribution  of  q 
by  qa(r,n,k),  one  may  invert  the  probability  inequality: 

P ^q0.05  (r»n’k)  < < q0.95  (r’n’k)]  = °-90  (26) 


to  obtain  a 905  confidence  interval  for  3 as: 

A A 

8/ qQ  95  (r,n,k)  < 3 < 8/qQ  05  (r,n,k) 

(27) 

A median  unbiased  estimate  of  g would  be: 

1 

A /\ 

3 — 3/ qQ  gy  ( > n > k ) 

(28) 

0 

11 

1 

1 i 

| 

--  , T --  --  - r — — — - 

J 

— 


As  proposed  in  [6]  the  ratio  R of  the  upper  to  lower  ends  of  a 
confidence  interval  on  the  shape  parameter  is  a useful  measure 
of  the  precision  with  which  the  shape  parameter  is  determined 
in  a sample  of  given  size. 

Arbitrarily  using  a 90%  interval  it  follows  from  Eq.  (27)  that 

R = q0.95/q0.05  (29) 

For  setting  confidence  intervals  on  Y it  is  noted  that  since 
the  random  variables  s and  q are  distributed  independently  of 
the  Weibull  parameters,  so  is  their  product  w*  = s • q,  i.e.: 

w*  = (y-y)3  *3/3  = (Y'Y)S  (30) 

Thus  from  the  distribution  of  w*(r,n,k)  one  may  set  a 90% 
confidence  interval  for  Y as: 

* wg  95  (r,n,k)  ^ w?  q (r,n,k) 

Y - — : * < Y < Y - : * (31) 

3 3 


A median  unbiased 

/\  A 

Y>  = Y 


estimate  of  Y is 

w0 . 50  (r’n’k) 

- 


3 


calculable 


as  : 


(32) 


A convenient  measure  of  the  precision  of  determination  for 
Y is  the  median  length  LQ  5Q  of  a (say)  90%  confidence  interval. 


From  Eq. 
as : 


(35)  and  the  definition  of  q,  Lq  may  be  calculated 


w*  - w* 
w0 . 95  w0 . 05 


J0 . 50 


3 q 


0.50 


(33) 
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Lq  50  depends  on  the  sample  size  parameters,  i.e.  k,  n and  r 
and  also  on  the  true  but  unknown  value  of  the  shape  parameter 
6. 


The  ratio  of  LQ  values  for  various  sample  size  choices  is 
however  independent  of  3. 


Substituting  Eq . (23)  into  Eq . (6)  and  using  Eqs.  (11)  and 
(12)  gives  the  following  expression  for  xp^  in  terms  of  the 


Pi 


= S 


;Y  4 -iogu  }8/b) 

i=l  3 = 1 


1/3 


(34) 


where  kp  = log(j4p) 


(35) 


The  population  value  of  xp^  is: 

c . 

Pi 


x . = k 1/6  s:Yn 

“ - n -i 


l o 


(36) 


^Dividing  Eq.  (34)  by  Eq.  (36)  and  raising  both  sides  to  the 
3-th  power  gives : 


x ~ ]1  *3/3  ~ * v n ^ ^ 


(37) 


The  right  hand  side  of  Eq.  (37)  involves  only  the  u^  and  the 
pivotal  functions  q and  w*  and  thus  the  function  on  the  left 
side  of  Eq.  (37)  or  its  logarithm  is  a pivotal  function. 
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Defining 


u*  (r,n,k,p)  = 3 log  [5^—]  (38) 

Pi 

one  can,  given  the  percentage  points  of  u* , set  confidence 

limits  on  x . . Two  sided  90%  confidence  limits  would  have  the 
Pi 

form : 


V • «p  [-“0.95/e1  ‘ xPi  •=  xPi  ' exp  ['uo.os/e]  (39) 

A median  unbiased  estimate  of  x . may  be  calculated  as: 

pi 

xpi  ■ xp  ‘ exp  [‘u0.50/61  140:1 

It  should  be  noted  that  all  of  the  foregoing  applies  even  when 
the  stress  at  which  it  is  desired  to  estimate  Xp  is  not  one 
of  the  stresses  at  which  life  tests  are  run.  Confidence  limits 
and  median  unbiased  point  estimates  can  be  calculated  for  any 
stress  within  or  without  the  range  of  stresses  encompassed  by 
life  tests. 

A measure  proposed  in  [6]  of  the  precision  with  which  Xp^  is 
determined  is  the  median  ratio  RQ  of  the  upper  to  lower 
ends  of  its  confidence  interval.  Choosing  a 90%  interval  this 
ratio  becomes,  from  (39) 


Ro.so  = exp  [- 


"U0. 05  + U0 . 95 


*0.50 


wherein  the  median  value  of  3 is  expressed  as  the  product  8q0  50 
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5.  NUMERICAL  SOLUTION  OF  THE  ML  EQUATIONS 


The  solution  of  Eqs.  (9)  and  (7)  by  means  of  a general  nonlinear 
equation  solver  routine  that  computes  derivatives  using  finite 
differences,  was  found  to  be  slow  and  unreliably  convergent. 


The  ad  hoc  solution  scheme  described  below  was  found  to  be 
quite  fast  and  dependable. 

Using  the  transformation 


yij 5 si  ' xi«) 


(42) 


Eq.  (9)  reduces  to 


i/e 


k ri 


i-  n • 

k i 


k ni 


+ ^ £ logy  / E r.  - E E y*  logy  / E E y?. 

i-1  j-1  1J  i=l  1 i=l  j-1  1J  1J  i=l  j-1  1J 


(43) 


k ni 


a 

+ y E E y?.  logS./  E E yj.  -y  E r.logS ./* 
i = l j-1  1J  1 i = l j-1  1J  1 x.“Ji 


= 0 


while  Eq.  (7) 


k 

E r . logS.  - 
i=l  1 1 


becomes 


n . 


E r - E E y . . , c 
i-1  1 i-1  j-1  1J  l0gSi 


k 

E 

i=l 


n . 

l 

E 

j-1 


1 

A 

yij 


0 


(44) 


When  Eq.  (44)  is  satisfied  the  last  two 
also  satisfied.  The  remaining  terms  of 
as  the  terms  that  must  sum  to  zero  when 
shape  parameter  from  k type  II  censored 


terms  of  Eq.  (43)  are 
Eq.  (42)  are  recognized 
estimating  the  Weibull 
data  groups  with  y^ 
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representing  the  j-th  ordered  observation  in  the  i-th  group 

(cf.  [7]).  For  a given  set  of  y . . the  solution  of  these 
* 1 J 

equations  for  3 is  readily  accomplished  by  the  Newton-Raphson 
method . 

From  this  observation  the  following  approach  to  the  simultaneous 
solution  of  Eqs.  (48)  and  (49)  suggests  itself: 

A A 

1)  Guess  a value  of  Y = yloi 

2)  Transform  the^data  to  y^  using  Y = Y in  Eq.  (47) 

3)  Solve  for  B =3  ^ from  the  first  three  terms  of  Eq.  (48) 

using  y.  . = y.1?^ 

4)  Use^B(0,and  Y(0'in  Eq.  (49)  and  solve  for  an  improved 
Y =Y  ^ using  the  Newton-Raphson  method* 

5)  Repeat  steps  (l)-(4),  iteratively  replacing  Y by  the 

latest  value  emerging  from  step  (4) , until  successive 

/\  -A 

Y and  3 estimates  do  not  differ  by  more  than  a 
prescribed  amount. 


This  procedure  was  incorporated  into  a simulation  program 
written  in  Fortran  IV,  that  generates  10  000  realizations  of 
a random  experiment  in  which  k samples  of  size  n are  tested  at 
each  of  k stresses  , S2 . • . until  the  first  r failures 
occur. 


The  failure  lives  follow  a two  parameter  Weibull  distribution 


with  shape  parameter  3 = 1.0  and  with  a p-th  percentile  value 


Xpi  at  stress  S^(i=l...k). 

The  input  data  consists  of  the 
and  x • (i=l,k)  and  the  values 

pi 

Sk+1’  Sk+2 ' ' ' Sm  at  which  U is 


values  n,  r,  k,S..  (i  = l,k),  p 
of  three  additional  stresses 
of  interest  to  calculate  Xpj- 
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By  setting  xpi  constant  for  all  i one  simulates  the  case  where 
the  stress-life  exponent  has  the  value  y = 0. 

A A A 

The  program  calculates  y,  $ and  xpi  (i -1 , . . . k+m)  for  each 
realization  and  compiles,  for  specified  functions  of  these 
observations,  their  empirical  distribution  over  10  000  simulated 
realizations  of  the  experiment. 

6.  NUMERICAL  EXAMPLE 

The  following  data  are  the  ordered  times  to  failure  in  rolling 
contact  fatigue  of  ten  hardened  steel  specimens  tested  at 
each  of  four  values  of  the  contact  stress. 


Stress 

10^-psi 

Ordered 

Lives 

0.87 

1.67,  2 

.20,  2.51,  3.00,  3.90,  4.70,  7.53, 

14.70, 

27.76  , 

37.4 

0.99 

0.80,  1 

.00,  1.37,  2.25,  2.95,  3.70,  6.07, 

6.65, 

7.05, 

7.37 

1.09 

0.012  , 

0.18,  0.20,  0.24,  0.26,  0.32,  0.32 

, 0.42, 

0.44, 

0.88 

1.18 

0.073, 

0.098,  0.117,  0.135,  0.175,  0.262, 

0.270, 

0.350, 

0.386,  0.456 

Rolling  contact  fatigue  data  of  this  type  are  customarily  treate 
as  samples  from  two-parameter  Weibull  populations  having  a 
scale  parameter  that  varies  inversely  with  a power  of  the 
contact  stress  and  a shape  parameter  that  is  invariant  with 
stress  [cf.  Lieblein  and  Zelen  (8)]. 
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For  the  most  part  these  assumptions  appear  to  be  satisfied  by 
the  present  data.  Probability  plots  suggest  however  that  the 
first  failure  at  the  stress  S = 1.09  is  an  outlier  and  that 
at  S = 0.87  a three  parameter  Weibull  model  with  a location 
parameter  of  about  1.50  is  indicated. 

For  expository  purposes  we  nonetheless  accept  the  two -parameter 
Weibull  model  as  adequately  representing  the  data. 

To  test  the  constant  shape  parameter  assumption  we  calculate 
the  raw  ML  shape  parameter  estimates  for  each  sample,  obtaining 
the  values  listed  below: 

Stress  (106  psi)  0.87  0.99  1.09  1.18 

ML  Shape  Parameter  Est.  0.953  1.57  1.43  1.96 

Following  [9]  we  employ  as  test  statistic  the  ratio  1.96/0.953  = 

2.06  of  the  extremal  shape  parameter  estimates. 

From  values  tabled  in  [9]  we  find  that  for  k = 4 , n = r = 10, 
the  90-th  percentile  of  the  null  distribution  of  the  extremal 

/\  /v 

shape  parameter  ratio  is  ,, ) n Qn  = 2.47.  Since  2.06<2.47 

iTicix  mill  u.yu 

there  is  no  reason  to  reject  the  hypothesis  of  a common  shape 
parameter. 

Proceeding  then,  we  calculate  the  joint  ML  estimates  of  the 
stress -life  exponent  y and  the  common  shape  parameter  3 from 
Eqs.  (8)  and  (10)  to  be: 

n 1 

Y = 13.89 
3 = 1.166 
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Additionally  we  calculate  the  ML  estimate  3 of  the  common 
shape  parameter  3 that  applies  when  the  are  not  constrained 
in  any  way  (cf.  [7])  with  the  result  3 (1)  = 1.343. 

A 

The  estimate  3 is  closer  to  the  average  of  the  four  individual 
ML  estimates  tabled  above  than  is  the  estimate  3 = 1.166  that 
results  when  the  scale  parameters  are  constrained  to  vary  with 
an  inverse  power  of  stress.  This  suggests  some  lack-of-fit  to 
the  inverse  power  law. 

The  ML  tenth  percentile  estimates  at  the  four  stresses  are 
as  follows : 

Stress  CIO6  psi)  0.87  0.99  1.09  1.18 

~ power  law  2.209  0.3672  0.0965  0.0321 

unconstrained  2.328  0.7862  0.0656  0.0457 

The  simulation  program  was  run  with  n = r = 10,  k=4  and  = 
0.87,  S2  = 0.99,  = 1.09,  S^  = 1.18  and  the  additional  value 

S^  = 0.75.  The  distributions  were  calculated  for  q,  w* , and 
the  u*  applicable  at  each  S^.  The  5-th,  50-th  and 
95-th  percentiles  of  these  random  variables  are  listed  in 
Table  1. 

Using  these  tabled  values  one  calculates  from  Eq . (27)  the  90% 
confidence  interval 

0.913  = 1.166/1.277  < 3 <1.166/0.8459  = 1.378 

From  Eq.  (28)  a median  unbiased  estimate  of  the  shape  parameter 
is 

A 

3’  = 1.166/1.024  = 1.139 
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From  Eq.  (29)  the  precision  measure  R is  calculated  as 

= 1.277  = 

K (5T8459 

This  value  is  in  good  agreement  with  the  value  calculated  for  a 
single  uncensored  sample  of  size  n = 40,  indicating  that  there 
is  a negligible  loss  in  precision  for  estimating  8 by  conducting 
tests  at  four  stresses  rather  than  at  only  one. 


From  Eq.  (31)  a 90%  confidence  interval  for  the  stress-life 
parameter  is  calculated  to  be 


11.92  = 13.889  - 2.293/1.166  <Y<  13.889  + 


2.433/1.166  = 15.98 


A median  unbiased  estimate  of  y is  from  Eq . (32) 


Y'  = 13.889  + .3783/1.166  = 14.21 


From  Eq.  (33)  the  product  of  8 and  the  precision  measure  Lq 


'0.50 


2.293  + 2.433 

T7V21 


4.62 


Using  Eqs . (39),  (40)  and  (41)  yields  the  following  values  of 
the  median  unbiased  estimates  and  90%  confidence  limits  for 
X«  m-,  i = 1....5.  Also  listed  for  each  stress  are  the  values 

u . 1 1)  1 g 


of  the  precision  measure  Rq 
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Stress 

Median  Unbiased 

Xq  ^q  Estimate 

901  Confidence 
Lower  X°*10 

Interval 

Upper 

d3 

-0.50- 

0.75 

16.55 

7.22 

38.4 

6.70 

0.87 

2.13 

1.09 

3.86 

4.21 

0.99 

0.358 

0.193 

0.572 

3.42 

1.09 

0.094 

0.050 

0.152 

3.55 

1.18 

0.031 

0.016 

0.054 

4.03 

Figure  1 shows  on  logarithmic  scales  the  straight  line  that  joins 
the  raw  estimates  of  xQ  plotted  against  stress  along  with 
the  bands  formed  by  the  upper  and  lower  confidence  limits  cal- 
culated at  each  stress. 


Also  shown  are  the  .-iL  estimates  of  Xq  at.  the  four  stresses 
at  which  data  were  taken  computed  under  the  assumption  of  a 
common  shape  parameter  but  with  no  constraint  on  the  variation 
of  Xq  with  stress  and  their  associated  90%  confidence  inter- 
vals . 


The  lower  confidence  limit  for  the  uncontrained  estimate  at  = 
0.99  just  touches  the  line  fitted  under  the  power  law  constraint. 
This  confirms  the  indication  of  lack  of  fit  to  the  power  law 
model  suggested  by  the  comparison  of  the  constrained  and  un- 
constrained shape  parameter  estimates  discussed  earlier. 


Fig.  (2)  is  a plot  of  Rq  against  stress  level.  Fig.  (2) 
illustrates  how  Rq  ^.q  is  a minimum  near  the  midpoint  of  the 
stress  range,  increasing  substantially  when  extrapolating  to 
stresses  outside  the  range  of  the  tests. 
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The  value  of  5q  for  a single  test  of  size  n =r  = 40  was 
found  by  simulation  to  be  = 3.41.  This  compares  favorably 

with  the  minimum  value  shown  on  Fig.  2.  This  indicates  that  if 
high  precision  is  required  at  some  specific  stress  level  a 
negligible  loss  in  precision  is  suffered  if,  rather  than  con- 
ducting all  tests  at  this  stress  level,  some  specimens  are  tested 
at  surrounding  stress  levels.  Plots  like  Fig.  (2)  can  be 
constructed  in  advance  of  any  actual  testing  to  determine  the 
choice  of  stress  levels,  the  number  of  levels  and  the  sample 
sizes  that  yield  a precision  "profile"  that  the  experimenter  finds 
suitable.  In  addition  as  noted  previously  the  sample  size  and 
censoring  number  can  differ  at  each  stress  level. 


1 


ACKNOWLEDGEMENT 


Derivation  of  the  ML  estimation  equations  and  pivotal  functions 
was  performed  under  contract  to  Mobil  Research  and  Development 
Co.  under  the  project  monitorship  of  Dr.  P.  E.  Fowlcs.  This 
support  and  permission  to  publish  this  portion  of  the  work  is 
gratefully  acknowledged.  Development  of  a solution  scheme  and 
all  simulations  and  numerical  work  was  performed  under  the  spon- 
sorship of  the  Air  Force  Office  of  Scientific  Research  under 
Contract  No.  F49620 - 77 -C-0007  with  Dr.  I.  N.  Shimi  as  Project 
Monitor.  This  support  is  likewise  gratefully  acknowledged. 


TABLE  1 


q = 3/6 

A A 

w*  = (y-y)6 

A A 

u*  = Blog  [x 

S = 0.75 
S = = 0 

S - s2  - 0 

S = S,  = 1 

S - s43  - 1 


Percentiles  of  Pivotal  Functions 
k=4,n=r=10 

S.-0.87,  S2=0.99,  S3=1.09,  S4=1.18 


0.05 

0.50 

0.95 

0.8459 

1.024 

1.277 

-2.433 

-0.3783 

2.293 

0.10/x0.10^ 
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.87 
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0.0441 
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.99 
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.09 
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0.0318 
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.18 
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0.0396 
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